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ABSTRACT 

Real root finding of polynomial equations is a basic problem 
in computer algebra. This task is usually divided into two 
parts: isolation and refinement. In this paper, we propose 
two algorithms LZ1 and LZ2 to refine real roots of univariate 
polynomial equations. Our algorithms combine Newton's 
method and the secant method to bound the unique solution 
in an interval of a monotonic convex isolation (MCI) of a 
polynomial, and have quadratic and cubic convergence rates, 
respectively. To avoid the swell of coefficients and speed 
up the computation, we implement the two algorithms by 
using the floating-point interval method in Maplel5 with the 
package intpakX. Experiments show that our methods are 
effective and much faster than the function RefineBox in 
the software Maplel5 on benchmark polynomials. 



Keywords 

Real root refinement, Newton's method, floating-point, 
interval method 



1. INTRODUCTION 

Solving for real solutions of univariate polynomial equa- 
tions is a fundamental problem in computer algebra. Many 
other problems in mathematics or other fields can be re- 
duced to it such as real solving multivariate polynomial 
equations |23l El [3 [5], studying the topologies of real al- 
gebraic plane curves [Hj, and generating ray-traced images 
of implicit surfaces in computer graphics [17] ■ There is a 
vast literature on calculating real zeros of univariate poly- 
nomials. We refer to 20 for some of the references. 

The process of reliable computing real roots is usually di- 
vided into two steps: real root isolation (i.e., cutting the real 
axis so that each real root is contained in a separate inter- 
val) and real root refinement (i.e., narrowing each isolating 
interval to a given width). We mainly concerned with how 
to efficiently refine each real root of a univariate polynomial 
equation to a high precision. We propose two algorithms 
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LZ1 (Algorithm [3| and LZ2 (Algorithm 0} to refine real 
roots. 

LZ1 is a combination of Newton's method and the secant 
method. It is based on Theorem 12.11 (cf. Theorem 4.6 in 
[34]) that can help choose a starting point and guarantee the 
quadratic convergence of the point sequence for Newton's 
method in an interval. By this theorem, all of the points 
in the sequence locate on the same side of the real root £, 
i.e., each point in the sequence provides a bound of £ of the 
same kind (upper bound or lower bound). For each point 
in the sequence, to get a bound of £ of the other kind, the 
secant method should be applied. Theorem 15.21 shows that 
LZ1 has an at least quadratic convergence rate. 

LZ2 is also a combination of Newton's method and the 
secant method. But it makes an opposite choice of the start- 
ing point for Newton's method to LZ1. As a result, if we 
consider the starting point as a bound of £, then the point 
found by Newton's method becomes another bound of £ of 
the other kind. After that, by the secant method, a new 
bound of £ of the same kind with the starting point can be 
obtained. In this way, the two kinds of bounds of £ arise al- 
ternately. Thus, the precision of each bound of £ can benefit 
from its predecessor in this point sequence. So, it is not sur- 
prising that the sequence of inclusion intervals consisting of 
these bounds converges to £ at least cubically (cf. Theorem 

IQJi . 

However, the implementations of LZ1 and LZ2 will be- 
come slow if we evaluate polynomials in LZ1 and LZ2 ex- 
actly, since the representations of the numbers will swell dra- 
matically. Fortunately, the floating-point interval method 
can help solve this problem. But this method can only deal 
with the so called well-posed problems in numerical com- 
putation. In fact, to apply this method to speed up the 
computation, we have already treated all the ill-posed cases 
by exact methods before calling LZ1 or LZ2. 

To this aim, a square-free decomposition should be done 
first to get a list of square-free polynomials as components 
and to know the multiplicity of corresponding roots in each 
component. Then, we make a local monotonic convex de- 
composition (LMCD, cf. Definition [2} for each of these 
square-free polynomials to make sure each component has a 
monotonic convex isolation (MCI, cf. Definition [Q. After 
that, based on existing methods for real root isolation [9] [2] 
[32] EU El [HI [28] IU El], we compute a MCI for any polyno- 
mial in a LMCD of every square-free polynomial obtained 
in the first step. Finally, floating-point interval versions of 
LZ1 and LZ2 can be called safely in theory. 

We have implemented LZ1 and LZ2 in Maplel5 with the 



package intpakX which contains functions for floating-point 
interval computation in arbitrary precisions. Experiments 
were done on Chebyshev polynomials of the first kind to 
compare the efficiencies of LZ1, LZ2 and the Maple function 
Ref ineBox. The timings show that our implementations of 
LZ1 and LZ2 are much faster than Ref ineBox and that LZ2 
is usually faster than LZ1. 

Some variants of Newton's method also exist to reliably 
refine real roots of univariate polynomial equations. 

Ramon E. Moore [IS] in 1966 gave an interval Newton's 
method to compute a real root of a univariate real function 
in a closed interval. This method can be speeded up by using 
floating-point interval arithmetics [27]. It can give verified 
results, but sometimes it may fail, e.g. the resulting interval 
may contain the original one. 

George E. Collins and Werner Krandick 11 in 1993 pre- 
sented three versions of a variant of Newton's method and 
proved that the exact version has quadratic convergence 
rate. Among the three algorithms, they showed that the one 
using floating-point interval arithmetics performed much ef- 
ficiently than the other two. Note that the restrictions of 
the input inclusion interval and the selection of the start- 
ing point for Newton's method in the exact version of their 
method are also the same with that in Theorem 12.11 But 
instead of a secant, when it is needed in LZ1, they used 
the line parallel to the tangent at the point that Newton's 
method is applied. Moreover, their method needs the input 
interval ro satisfy a "3/4-assumption" to get started. 

John Abbott [T] in 2006 proposed an algorithm named 
QIR. Then, based on the QIR, Michael Kerber and Michael 
Sagraloff [15] in 2011 gave an algorithm EQIR and analyzed 
its bit complexity. QIR and EQIR are combinations of the 
secant method and the bisection method. One feature of 
them is that they do not need any information of derivatives. 

The rest of the paper is structured as follows. In Sec- 
tion [21 we introduce the basic theorem that LZ1 is based 
on. Section [3] is devoted to describing how multiplicities 
and the square-free polynomial components should be cal- 
culated. The notions of MCI and LMCD are given in Section 
[4] in order to get the intervals that can be used as inputs of 
LZ1 and LZ2. In Section [S] the main algorithms LZ1 and 
LZ2 are presented and their superlinear convergence prop- 
erties are proved. Then we discuss how to use floating-point 
interval method to speed up LZ1 and LZ2 in Section [5] At 
last we make comparisons on efficiency among LZ1, LZ2 and 
the Maple function Ref ineBox in Section [7] 

2. BASICS 

This section is about a theorem which provides the basis 
of Algorithm [3] named LZ1 in Section [5] 

Newton's method is a fast algorithm for computing solu- 
tions of equations. For a function and a stating point, gen- 
erally, it is not easy to determine whether Newton's method 
converges or not. However, for univariate real functions 
there exists such a theorem as follows (cf. Theorem 4.6 in 
[34] written by Qingjin Yan in 1992). 

Theorem 2.1. If a function f £ C 2 [a,b] and satisfies the 
following conditions: 

1. f(a)f(b) < 0; 

2. the sign of f"(x) does not vary on [a,b\; 



3. f'(x) / when x £ [a, b]; 

4. x G [a, b], f(x )f"{x ) > 0, 
then the sequence {xk} generated by 

x k+1 = x k - ; \ , k = 0, 1, . . . 

f'{Xk) 

monotonically converges to the unique real root of f(x) = 
in [a,b] and the convergence rate is at least 2. 

Remark 2.2. Note that the second condition of Theorem 
\2.1\ allows f"(x) to reach zero when x £ [a,b]. If we re- 
strict this condition to "f"{x) 7^ when x £ [a, b] ", then the 
iteration sequence will strictly monotonically converges. 

3. MULTIPLICITIES 

We briefly recall how to compute multiplicities of roots of 
polynomial equations by existing methods. 

Though numerical tools exist (e.g. [291 116] ) for comput- 
ing multiplicities of roots of polynomial equations and most 
probably their outputs are correct, it is still possible that 
they output incorrect results sometimes. In other words, for 
a univariate polynomial / £ Q[x], computing the multiplic- 
ities of real or complex roots of the equation / = is an 
ill-posed problem in numerical computation [27| . 

Thus, to know these multiplicities exactly, we should use 
exact methods. We can make a square-free decomposition of 
/ such that there exist s polynomials gi, . . . , g s in Q[x] with 
f = gl 1 . . . g r s s where 1 < ri < ■ ■ ■ < r s < n and (gi,gj) = 1 
when i 7^ j. Then f — has a (complex or real) root £ with 
multiplicity r if and only if there exists a polynomial gi in 
the square-free decomposition of / such that r = and £ 
is a single root of gi — 0. The square-free decomposition 
of a univariate polynomial is easy to compute in the soft- 
ware Maplel5 with the function sqrfree. A factorization 
with the function factor in the same software can also help 
find the multiplicities of corresponding polynomial zeros. In 
practice, a square-free decomposition or a factorization costs 
much less time than real root isolations and refinements and 
can help reduce an reducible polynomial in Q[x] with high 
degree to some lower degree polynomials which are relative 
easier to deal with. 

In the rest of this paper, we assume that / is a square-free 
polynomial with rational coefficients. 

4. MONOTONIC CONVEX ISOLATIONS 

To refine a real root of a square- free polynomial / G Q[a;], 
it should be first isolated from other real roots. In this sec- 
tion, we give a concept called monotonic convex isolation 
(Definition 14. ip . If / has such an isolation, then any real 
root of / is contained in a closed interval where the mono- 
tonicity and convexity of / is fixed unless the interval is a 
point. The refinements of real roots in the next section will 
benefit a lot from the properties of monotonicity and con- 
vexity of / on these isolating intervals. However, certain 
square- free polynomials in Q[x] have no monotonic convex 
isolations. In this case, we want to decompose / to a multi- 
plication of some polynomials that have monotonic convex 
isolations (Theorem I4.5[) . This leads to the concept of local 
monotonic convex decomposition (Definition 14. 4[) . At the 
end of this section, we give an algorithm to compute such a 
decomposition. 



Definition 4.1 (MCI). Given a real root isolation I so 
of a square-free polynomial f £ Q[x], we say that I so is a 
monotonia convex isolation of f if I £ Iso is not a point 
implies that f'(x) 7^ and f"(x) 7^ for all x £ I . 

Remark 4.2. "I £ Iso is not a point" in the above def- 
inition means that the closed interval I is not in the form 
of [a, a] where a is a rational number. Note that such stuff 
may exist in Iso. 

If / has a MCI, then we can work it out by Algorithm [T] 
Note that Algorithm [1] is not a definite algorithm but a de- 
scription of a class of algorithms. When the base algorithm 
for real root isolation (e.g. the RS C-library [221 1231 124| of 
Fabrice Rouillier and DISCOVERER [H [31] of Bican Xia) is 
chosen, then Algorithm [T] will become definite. Hence, its 
efficiency mainly depends on the chosen base algorithm. 

Algorithm 1: MCI 

Data: A polynomial / £ Q[x] that has a MCI 
Result: A monotonic convex isolation of / 

1 Compute a real root isolation Iso of / in common sense 
by a base algorithm. 

2 For each interval I £ Iso that is not a point, cut I by 
bisection method and consider the sign of the value of / 
at the midpoint of I to find the interval that contains 
the real root in I (sometimes, the midpoint is just the 
real root). 

3 Repeat step 2 until neither /' nor /" has real roots in 
the resulting closed interval I* (this can be done by 
using inclusion and exclusion criteria). 

4 Return the set consisting of all the 7* in step 3 and all 
the point intervals arising in Iso or during the 
computation of step 3. 

It is obvious that for a nonconstant square- free polynomial 
/, if it satisfies (/, /") = 1 then / has a MCI. However, not 
all real polynomials have monotonic convex isolations. We 
give a counterexample below. 

Example 4.3. Pick f as x{x + l)(x + 2). Then f = 
3x 2 + 6x + 2 and f" — 6(x + 1). It is easy to see that —1 
is the common real root of f and f" . Obviously, f has no 
MCI. 

Hence, we want a decomposition of / such that each com- 
ponent has a MCI. This idea delivers the following definition. 

Definition 4.4 (LMCD). // a square-free polynomial 
f £ Q[ie] can be represented as the multiplication of s poly- 
nomials gi, . . . , g B where deg(<?i) = 1 or gcd(g;, gr") = 1 for 
every i — 1, . . . , s, then we say that f has a local monotonic 
convex decomposition. 

Theorem 4.5. Every nonconstant square-free polynomial 
f £ Q[x] has a local monotonic convex decomposition. 

Proof. For a polynomial / £ Q[x] \ Q, we prove the 
theorem by induction on the degree d of /. 

When d — 1, the theorem holds obviously. Supposing 
the conclusion of the theorem holds when d < k (k > 1), we 
prove that it also holds for d = fc+1. Denote g := gcd(/, /"). 
It is easy to see that deg(g) < deg(/) when d > 1. If 
g = 1 then the conclusion holds. Otherwise, there exists 



a polynomial h £ Q[x] \ Q such that / = gh. Since deg(gr) 
and deg(h) are all no larger than k, we know that g and 
h all have local monotonic convex decompositions by the 
induction assumption. Thus, / has a local monotonic convex 
decomposition. 

Therefore, the conclusion of the theorem holds for all non- 
constant square-free polynomials in Q[x]. □ 

Corollary 4.6. Any nonconstant irreducible polynomial 
f G Q[x] forms a LMCD of itself. 

Remark 4.7. Sometimes, the irreducibility of a polyno- 
mial f in Q[x] is easy to test by certain irreducibility criteria, 
e.g. Etsenstein's criterion. If the irreducible decomposition 
of f i' n Q[ x ] * s known, then it is also a monotonic convex 
decomposition of f according to Corollary \4- 



Based on the proof of Theorem 14.51 we give Algorithm [2] 
(LMCD) to compute a local monotonic convex decomposi- 
tion of /. The termination and correctness are obvious and 
we omit the proofs. In most cases, the algorithm LMCD 
only tests that / and /" are coprime. This test is usually 
very fast in practice. 

Algorithm 2: LMCD 

Data: A nonconstant square-free polynomial / £ Q[x] 
Result: A local monotonic convex decomposition of / 

1 S:={}; 

2 if deg(/) = 1 then 

3 S :={/}; 

4 return S; 

5 end 

6 g:= gcd(/,/"); 

7 if g = 1 then 

8 S :={/}; 

9 return S; 

10 end 

11 S := S U LMCD( ff ) U LMCD(//g); 

12 return 5*; 



Example 4.8. For the polynomial f = x 3 + 3x 2 + 2x in 
Example \4-3\ it can be decomposed into two polynomials by 
Algorithm^ i.e., the output is {x+l,x 2 +2x}. Forx+l = 0, 
we can directly compute the root. For x 2 + 2x, it has a 
monotonic convex isolation. 



5. REAL ROOT REFINEMENTS 

In this section, we study how Newton's method and the 
secant method can be combined and applied on the MCI 
intervals to compute narrower inclusion intervals of the real 
roots of a univariate polynomial equation. For this, we pro- 
vide two algorithms LZ1 (Algorithm [3| and LZ2 (Algorithm 
[4]), and prove that their convergence rates are at leat 2 and 
3, respectively. 

We design LZ1 (Algorithm [3]) based on Theorem [2Tl Ac- 
cording to this theorem, the initial point xo for Newton's 
method should satisfy the condition f(xo)f"(xo) > 0, and 
the point sequence {xi}°l of Newton's method converges 
monotonically. This provides a series of bounds of the real 
root ^ on one side. Since the convexity does not change in 
the input MCI interval [a, b], the secant method can be used 



to obtain a bound d+i on the other side of £ after getting 
each Xi by Newton's method. Then a sequence of inclu- 
sion intervals of £ is obtained. Theorem l5.2l shows that this 
sequence converges at least quadratically. 

The second main algorithm LZ2 (Algorithm 2| is also a 
combination of Newton's method and the secant method. 
In contrast to LZ1, the condition for selecting initial point 
Co for Newton's iteration in LZ2 is /(co)/"(co) < 0. As a 
result, the point z {x\) obtained by Newton's method on / 
at Co locates at the other side of £. Then we can get c\ as a 
bound of £ at the same side with Co by the secant method. In 
this way, the inclusion interval sequence of £ can obtain an 
at least cubic convergence rate (cf. Theorem I5.4|l . However, 
note that when z £ [a, b], the behavior of / at z will be out 
of control. In this case, [a, 6] should be cut until the new z 
locates in the new initial MCI interval. This task is easy to 
do by using the secant method as described in lines 1 11 II of 
Algorithm [4] Then the fast convergent iteration can start. 

We want the outputs of LZ1 and LZ2 to be narrow in- 
tervals from which the precisions of the corresponding real 
roots can be read out directly. However, no matter how 
narrow an interval containing zero is, we cannot obtain any 
correct bits by comparing the two endpoints of the interval. 
Hence, LZ1 and LZ2 are not allowed to input intervals that 
contain zero. In fact, when a monotonic convex isolation 
of a square-free polynomial / is performed as we discussed 
in the last section, it is very easy to check whether f — 
has zero solution or not. If such solution exists, its isolation 
interval can be set as the interval [0, 0] in the MCI directly. 

Example 5.1. Given L = 8 (decimal digits), a square- 
free polynomial f := x 3 — 20a; + 7 and an interval [y§§§, ] 
in a monotonic convex isolation of f, we show the process 
of the iteration of Algorithm\3\ as follows. 



Algorithm 3: LZ1 



[ao, b 
[oi,6ij = [ 
[02,62] = [ 



4389 1097, 



1024 ' 256 ' 

40379863349 80788619485, 



9422150912 ' 18851042816 J 
18537520582609738516073933520140272183127 
4325505299448020965613362785214886187776 ' 
503844952756268930733017297728389, 



117566100791919345105185727965440 ' 

(62 - a 2 )/b 2 w 0.6053885328 x 10~ 14 < 10~ s . Hence, the 
algorithm terminates and outputs [02,62]. 

Theorem 5.2. The iteration in Algorithm [3| converges 
and the convergence rate is at least 2. 

PROOF. From the loop in lines 13131 of Algorithm [3] we 
know that the iteration is 



Xi+l = x 



Ci+l — 



f'(?i) 
XifjCj) - Cjf(Xi) 



(1) 

(2) 



f(*)-f(xi) 

where Xi and Ci means the i-th iteration of x and c for i = 
0, 1, . . ., respectively. 

We first prove that the real root £ of the equation / = 
always lies in the open interval between x% and d. By 
Theorem 12.11 we know that the sequence {xi}%jL converges 
monotonically to £ and that no Xi is equal to £. Thus, these 



Data: A closed interval [a, b] (ab > and a < 6) in a 
monotonic convex isolation of a square-free 
polynomial / £ Q[x] and a positive integer L 
Result: A closed interval [a',b'] such that the only real 
root £ of / = in [a, 6] is also in it and 
(6' - a')/|e| < W- L 



1 if \b-a\ < 10~^min(|c 

2 j return [a, 6]; 

3 end 

4 if f(a)f" (a) > then 

5 j x < — a; c < — 6; 

6 else 

7 j x < — 6; c < — a; 

8 end 

9 while 



then 



— c| > 10 L min(|a;|, |c|) do 
fix)', 

fi&, 

x - u/f'(x); 
(xv — cu)/{v — u) ; 

p; 

15 end 

16 return [min(a;, c), max(i, c)]. 



10 
11 
12 

13 
14 



u 

V ■ 



C i 
X 



Xi's are the bounds of £ on the same side. We only need 
to prove that the Ci's are the bounds of £ on the other side 
and d 7^ £. Note that if a and Xi are on different sides of 
£, then f {*)/(/(*) - f(xi)) and -f(xi)/(J(a) - /(*<)) are 
all in (0, 1) and their sum is 1. By induction on i, applying 
Jessen's inequality we have that /(ci+i) < when f"(a) > 
and /(ci+i) > when /"(a) < 0, i.e., 

/(a+i)/"(a) < 0; 

consequently, Cj+i locates on the opposite side of Xi+i w.r.t. 
£ and Ci+i 7^ £. This proves that £ always lies in the open 
interval between Xi and c;. 

Next, we study the convergence of the sequence consisting 
of the interval lengths. According to the Mean Value The- 
orem, there exist a A € (min{xi, £}, max{ij, £}) and a /1 G 
(min{ci,a;i},max{ci,a;i}) such that f(xi) = f'(X)(xi — £) 

and f(c ' c ) Z f x < ' x ' ) = /'(m). Note that the sequence {c»}g 
converges to £ monotonically. Indeed, since Cj+i and Ci are 
on the same side of £ and 



Ci+l 



a) fia) 



fid) - fiXi) 



fir) 



ia-0, 



we know that {ci}°% is monotone and has a bound £. Then 
{ci}^i converges and it is easy to see the limit is £. Since 
there exists ip between a and Xi such that f(ci) = f{xi) + 
(a - Xi)f'(xi) + (a - Xi) 2 f"(ip)/2, we have 



Xi+l 



Ci+l 



fiXi)( 



fjxi) _ XjfjCj) - CjfjXj) 
fiXi) fiCi)-f{Xi) 

1 1 , 



fixi)f"iv) 



2f'( t i)f'(x i ) 

f"{<p)f'W 



fiXi) 
(Xi — d) 

ixi - Ci)(xi - £). 



From Definition 14. II and the input of Algorithm [3] we know 
that /'(£) / and /"(£) / 0. Hence, 

i im supi^-^i <|il| /0 . 

Therefore, the sequence of the lengths of the intervals has 
at least quadratic convergence rate. □ 



Algorithm 4: LZ2 

Data: A closed interval [a, b] (ab > and a < b) in a 
monotonic convex isolation of a non-constant 
real univariate square-free polynomial / and a 
positive integer L 

Result: A closed interval [a', b'] such that the only 
root £ of / = in [a, b] is also in it and 

0' -a'VICI < 10~ L 

1 if |6 - o| < 1CT L min(jaj, |6|) then 

2 | return [a, b]; 

3 end 

4 if /(a)/"(a) > then 



5 




x < 


a; c •< — b; 


6 


e 


lse 




7 




X i 


— b; c i — a; 


8 


end 




9 


U i 


f(x); v < — /(c); 


10 


Z 




c-v/f'(c); 


11 


while 


z ^ [a, 6] do 


12 




C -t 


(l!) — Cu)/{v — 


13 




V i 


— /(c); 


14 




Z i 


— c-v/f(c); 



15 end 

16 x i — z; 

17 while \x — c\ > 10 

18 
19 
20 



— L . 



min(|x|, |cj) do 

u i — f(x); 

c < (xv — cu)/(v — u); 

if |x — c| < 10~ L min(|a;|, |c|) then 

21 | return [min(j;, c), max(i, c)]; 

22 end 

23 V < /(c); 

24 x < — c—v/f'(c); 

25 end 

26 return [min(x, c), max(x, c)]. 



Theorem 5.4. The iteration in Algorithm [^] converges 
and the convergence rate is at least 3. 

Proof. The iteration in Algorithm!?] has two stages, i.e., 
the loops in lines 14141 and in lines 14141 

We first prove that the loop in lines 3E] terminates. In 
this process, x is fixed. By Jessen's inequality, we have that 
f(d)f"(a) < for the i-th iteration of c in this loop (cf. the 
proof of Theorem I5.2|l . Hence, these Cj's all locate at the 
other side of x w.r.t. £ and if the loop does not terminate 
the sequence consisting of them monotonically converges to 
£. Thus, the sequence {zi}°Z will also converge to £ £ (a, 6), 
a contradiction. Therefore, the loop in lines 14141 terminates 
and the number of iterations is independent of L. 

Next, we study the loop in lines 14141 which is the iteration 

as 



Xi+l 



Ci+l = 



c m 

1 f'(a) 

Cif{Xi+l) - X i+1 f(Ci) 



f(d) 



(3) 
(4) 



where Co is the last value of c in the loop in lines |4]|4] and 
iro is the x assigned in lines 14141 Consequently, xi locates in 
the open interval between xo and £. By Jessen's inequality 
and from iteration (j4| we know that ci locates in the open 
interval between Co and £. Note that the derivative of the 
function x(t) :=*-/(*)//'(*) is f(t)f"(t)/(f (t)) 2 and that 
/(co)/"(co) < 0. By induction on i, it is not difficult to 
see that the two sequences {xi}°Z and {ci}^i au mono- 
tonically converge to £ but from opposite directions. This 
proves that the iteration in Algorithm [4] converges to £. 

At last, we study the convergence rate of this iteration. It 
is in fact the convergence rate of the second stage, i.e., the 
loop in lines [4]|4] Then we have that 



Xi+l — Ci+l 



f(Ci) af{x i+ i) - x i+1 f(a) 



f'(a) 



/(«)( 



/'(«)' 



Example 5.3. For the polynomial and the initial interval 
in Example \5.1l we apply Algorithm [4] to them and get the 
intermediate interval sequence (other than Examvle \5.1\ the 
intervals renew only one endpoint each time): 



[oo, bo] 
[ao,h] 
[oi, b{\ 
[ai, 62] 



r 4389 1097, 



1 1024 ' 256 ' 
r 4389 1261419417, 



1 1024 ' 294336896 J 

, 6671209230324943307293 1261419417 , 



1 1556645655550311117184 ' 294336896 
r 6671209230324943307293 



1 1556645655550311117184' 
283700456965465533230109 • • 



42 digits are omited 



66198056239039164770905 ■ • ■ 42 digits are omited 



(62 -ai)/b 2 » 0.1438320660 x 10~ 10 < 10" 8 . Hence, the 
algorithm Algorithm^terminates and outputs [01,62]- 



Expand f(xi+i) at a as f(ci)+f{ci)(xi+i-Ci)+f"(X)(xi+i- 
Ci) 2 /2 where n is in the open interval between c; and Xi+i. 
Then, 



Ci+l 



-fid) 



2/'(r?)/'( Ci ) 



f'(rff"(\) 
2/'(ij)/'(ci) 



where r £ (min(ci,^),max(ci,^)) and 77 G (min^+i, Cj), 
max(i, + i, Ci)). Hence, — Ci|}°^ and {\a — ^\}°1 Q have 
the same convergence rate around £. Now, we study the 
latter. Again, f(xi+i) can be expanded at c; as /(cj) + 
/ , (c i )(x i+ i-c i )+/"(c i )(x i+ i-c i ) 2 /2+/'"( K )( a;i+ i-c i ) 3 /6 
and /(^) (= 0) can be expanded as f(ci) + f'(ci)(^ — Cj) + 
/»(f-Ci) a /2 and /(c i )+/'(c0(e-c i ) + /"(ci)(e-c i ) 2 /2+ 



/"'(0)(£ - c 4 ) 3 /6. Then, from iteration Q we have that 
(c» - £)/(a!i+i) - - C)/(cO 



Ci+i - € = 



f(xi+i) - f{a) 
if (fit) + - «) + ^^(^+1 - 

=-pL(-/( c o+(^-$)/'(ci)- 
m f(Ci)(c ;)l /"'w/'w 2 f: a*) 

2/'(c0 /l K ^ 6/'( Cl ) 2 ( l V } 

1 , /"(<*)/» 
/'(r?) 1 4/'(ci) 

l,.,,,,., /'"(«)/'(r) 2 



6 

Therefore, 



~(.f"(0) + 



/'(Ci) 



))(«-€)* 



Ici+i-el _,3/"(0 2 + 4/' (£)/"' (0 



- lim I CI 3 



12/' (C) 2 



which implies that the sequence {d — £}£o nas at least cubic 
convergence rate, and so does the interval iteration using ([3]) 
and Q. □ 

Remark 5.5. The efficiency index of algorithms LZ1 and 
LZ2 are \/2 and y^3, respectively. Hence, in the sense of 
efficiency index LZ2 is also more efficient than LZ1. 

6. SPEED-UP 

In this section, we study how to speed up LZ1 and LZ2 
by using the floating-point interval method. 

As we have seen in Example 15.31 the sizes of the inte- 
gers representing the intervals increase dramatically. Given 
/ G Q[x] and p/q G Q, then f(p/q) G Q. In general, the 
maximal length of the numerator and denominator of f(p/q) 
approximates to deg(/) max(length(p), length(q)). Hence, if 
deg(/) is large, the computation will become quite time and 
memory consuming. 

Numerical computation with floating-point numbers can 
avoid this difficulty; however, we cannot know exactly whether 
a numerical result in common sense is reliable or not, i.e., 
its correct bits is unclear or whether it is zero is unclear. 

If floating-point numbers have definite representations, it 
is possible to estimate the bounds of a numerical result |26l 
1271 1141 IIP] , Hence, for a nonzero real number, a narrow 
interval with floating-point endpoints can show its correct 
bits. But zero recognition remains an ill-posed problem. 

Note that we have removed all the ill-posed cases in the 
former sections, i.e., the square-free decomposition to get 
multiplicities and remove common zeros of / anf /' (cf. Sec- 
tion^, the local monotonic convex decomposition to remove 
common zeros of / and /" (cf. Section [4]), the restriction 
ab > of input intervals of LZ1 and LZ2 to remove zero 
roots of / = (cf. Section[5|. As a result, the floating-point 
interval method can be applied to evaluate the bounds of 



f(p/q) (with Horner scheme) and decide its sign in LZ1 and 
LZ2. 

There exist several packages for the floating-point inter- 
val computation, e.g. Rump's INTLAB/Matlab package |27l 
125] . Geulig, Kraemer and Grimmer's intpack/Maple pack- 
age [13] based on the package intpak/Maple [12] developed 
by Corless and Connell, and Revol and Rouillier's MPFI 
open-source C library [22] , We use the functions in the 
intpack/Maple package, because this package can deal with 
floating-point numbers with arbitrary lengths and our algo- 
rithms need other functions in Maple. 

To apply the functions in intpacX, we should transform 
LZ1 and LZ2 into their floating-point interval versions and 
check the computation step by step. Since many similar 
details should be concerned, it is not suitable to present the 
pseudo codes here. So we give some principles as follows. 

1. Pick the initial length I (decimal digits) of floating- 
point numbers, e.g. take I := max(min(100, deg(/) + 
5), digits of a, digits of b). 

2. Replace every arithmetic in algorithms LZ1 and LZ2 
by corresponding interval arithmetic in intpakX. 

3. If an interval value of a polynomial at a point contains 
zero, then increase / (e.g. I := 21 for LZ1) and repeat 
relative computation until the endpoints of the new 
interval value share a fixed sign. 

4. Rewrite the inequality conditions in the loops and the 
first lines of the algorithms into interval versions, so 
that the algorithms can obtain stronger results than 
in the exact case. 

A natural question is whether the super-linear conver- 
gence rates can be retained or not when we use the interval 
method with floating-point numbers instead of exact com- 
putation with rational numbers. The answer is "YES" . The 
proof of this claim should deal with many details including 
why LZ1 and LZ2 can avoid ill-posed cases. We omit most 
of them and only show the key insight. Suppose that [u, v] 
is an interval containing £ and u, v are included in two nar- 
row enough intervals [u», it*] and [v*, v*], respectively, where 
it*, u*, w* and v* are all floating-point numbers. Then, 
^ G [«*,«*] C [it, v]. To determine whether the two intervals 
are enough narrow, we should check the signs of the values 
of the polynomial at the endpoints of [«,,«*] and [«*,«*]. 
However, [«*,«*] is a better choice than [«*,«»] in practice, 
since the former is less possible to obtain a false isolating 
interval and hence makes the algorithms more efficient than 
the latter. 

Example 6.1. Consider the numerical version of Exam- 
vle \5.3l Take the initial length of floats I := max(3 + 5, [1.6 x 
8]) = 12. Then, the numerical LZ2 yields 

[a ,b ] = [4.28515625000, 4.28613281250], 
[a ,bi] = [4.285 15625000, 4.285 631309351, 
[aiA] = [4.285631 226694662183, 4.285631 30935], 
[ai,b 2 ] = [4.285631226 694662183, 

4.285631226 709011277936552662771, and 
(62 - ox)/ai = 0.334818704119335811719986541959 x 10" 11 
< 10" 8 . 



To see the convergence rate, we continue to compute 

[02,62] = [4.285631226709011277936 4772441617276999 
1330501894, 4.28563122670901 1277936 55266277] . 

The correct digits of [a\, 61] and [02, 62] ar e 7 and 22, respec- 
tively. This shows that the iteration has cubic convergence 
rate. 

We can find that the expression of 62 in Example 16.11 is 
much shorter than that in Example 15.31 and that the swell 
of coefficients of polynomials has been solved. 

7. EXPERIMENTS 

To show the effectiveness of the algorithms LZ1 and LZ2, 
we implemented them in Maple 15 with the open source pack- 
age intpakX [13] and compared our implementations with 
the Maple function Ref ineBox in the package Regular/Chains 
[5] . All the experiments in this section were done on a com- 
puter with Intel(R) Core(TM) i3-2100 CPU @ 3.10GHz. 

Chebyshev polynomials of the first kind were used as the 
tested polynomials and were generated by the Maple func- 
tion ChebyshevT. The testing results are listed in Tables[T]2] 
In each table, L is the number of correct digits of the out- 
put (cf. Algorithms and Algorithm 2|, "ratiol" and "ra- 
tio2" are the rounded time ratios RefineBox/LZl and Re- 
fineBox/LZ2, respectively. For each polynomial, we aimed 
to refining the isolating interval to a precision 10 _i . The 
CPU times in the tables were tested by the Maple function 
time. 

For a polynomial / with degree n in Table [1] Table [3] 
and Table (4[ we isolated all its real zeros by the function 
RealRootlsolata^ 5] with the option ; rerr'=l/2, picked 
the (n/2)-th "box", and refined the isolating interval until 
its width was no larger than 10 -5 by Ref ineBosQ an d viewed 
the result as the initial isolating interval. For Table [3] and 
Tabled this result is [242345/262144, 484695/524288]. Each 
(n/2)-th "box" happened to contain a zero of the largest 
irreducible factor g in degree of / and contain no zeros of g' 
or g" . We denote the degree of this factor by n* . 

Factorizations for the polynomials in Table [2] were time- 
consuming. By Algorithm [21 we knew that these polyno- 
mials themselves compose their local monotonic convex de- 
compositions. Consequently, they all have monotonic con- 
vex isolations according to Definition 14.41 For a polynomial 
/ in Tabled we used Isolate (developed by Fabrice Rouil- 
lier in C language) in the RootFinding package instead of 
RealRootlsolate, because the former is faster than the lat- 
ter when n is large and Ref ineBox need not be called for the 
experiments in this table. The output intervals of Isolate 
with the option output=interval were narrow enough that 
/' and /" have no zeros and can be used as the input inter- 
vals of LZ1 and LZ2. We picked the (n/2)-th element in the 
output of Isolate. 

From all of these tables, we can see that our implementa- 
tions of LZ1 and LZ2 are much faster than the Maple func- 
tion Ref ineBox and in most cases LZ2 was more efficient 



There arc three other methods for isolating real roots in 
Maplcl5. They are realroot, KealRootlsolate with the option 
method= 'Discoverer' [32] [33] [30], and Isolate [22] [23] [24] in the 
RootFinding package. 

^According to Algorithm 6 in [5], this function performs a general- 
ization of the bisection method; but in this step, it took less than 1 
second in all tested cases. 



than LZ1. Since the environment variable Digits in our im- 
plementations of LZ1 and LZ2 did not change continuously, 
the time costs would have jumps in the tables. Hence, it 
is not confusing that sometimes LZ1 behaved a little better 
than LZ2 in our experiments. 



Table 1: Timings (s) (L = 1000) 



n 


n* 


RcfineBox 


LZ1 


ratiol 


LZ2 


ratio2 


100 


80 


89.560 


1.965 


46 


1.794 


50 


200 


160 


347.695 


4.492 


77 


3.900 


89 


300 


160 


349.317 


4.383 


80 


2.262 


154 


400 


320 


1414.819 


5.428 


261 


5.974 


237 


500 


400 


2232.982 


11.076 


202 


6.583 


309 


600 


320 


1427.814 


19.890 


72 


5.007 


285 


700 


480 


3370.198 


16.770 


201 


8.018 


420 


800 


640 


5898.132 


21.387 


276 


10.920 


540 


900 


480 


3332.259 


14.757 


226 


7.987 


417 


Table 2: More comparisons {L = 1000) 


n 


1000 


1200 


1400 


1600 


1800 


2000 


LZ1 


31.075 


27.284 


44.928 


27.346 


59.061 


34.398 


LZ2 


10.888 


18.626 


15.880 


25.958 


29.764 


33.275 



Table 3: Timings (s) (n = 1000, n* = 800) 



L 


RefineBox 


LZ1 


ratiol 


LZ2 


ratio2 


100 


55.629 


5.506 


10 


3.978 


14 


200 


217.075 


9.874 


22 


5.428 


40 


300 


540.262 


11.887 


45 


5.974 


90 


400 


1067.062 


11.934 


89 


5.787 


184 


500 


1809.892 


14.040 


129 


8.361 


216 


600 


2792.495 


13.915 


201 


12.480 


224 


700 


4103.359 


13.884 


296 


12.558 


327 



Table 4: More comparisons (n = 1000, n* = 800) 



L 


800 


900 


1000 


2000 


3000 


LZ1 


14.008 


19.578 


48.297 


48.578 


63.679 


LZ2 


12.448 


12.542 


12.542 


27.705 


27.612 



8. CONCLUSION 

In this paper, we provided a quadratically convergent al- 
gorithm LZ1 and a cubically convergent algorithm LZ2 to 
refine real roots of univariate polynomial equations. Before 
applying LZ1 or LZ2, the polynomial should be make the 
square-free decomposition and the local monotonic convex 
decomposition (LMCD), so that every component polyno- 
mial has a monotonic convex isolation (MCI). Moreover, we 
used the interval method with floating-point numbers to es- 
timate the values of polynomials at rational points to im- 
prove the efficiency the algorithms. Experiments on bench- 
mark polynomials showed that if we need high precisions of 
the real roots, then both of the two algorithms are much 
faster than the function RefineBox in Maplel5. 
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